/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.analysis.solvers;

import bd.org.apache.commons.math.exception.NoBracketingException;
import bd.org.apache.commons.math.util.FastMath;
import bd.org.apache.commons.math.util.Precision;

/**
 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
 * Brent algorithm</a> for finding zeros of real univariate functions.
 * The function should be continuous but not necessarily smooth.
 * The {@code solve} method returns a zero {@code x} of the function {@code f}
 * in the given interval {@code [a, b]} to within a tolerance
 * {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and
 * {@code t} is the absolute accuracy.
 * The given interval must bracket the root.
 *
 * @version $Id: BrentSolver.java 1244107 2012-02-14 16:17:55Z erans $
 */
public class BrentSolver extends AbstractUnivariateSolver {

    /**
     * Default absolute accuracy.
     */
    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;

    /**
     * Construct a solver with default accuracy (1e-6).
     */
    public BrentSolver() {

        this(DEFAULT_ABSOLUTE_ACCURACY);
    }

    /**
     * Construct a solver.
     *
     * @param absoluteAccuracy Absolute accuracy.
     */
    public BrentSolver(double absoluteAccuracy) {

        super(absoluteAccuracy);
    }

    /**
     * Construct a solver.
     *
     * @param relativeAccuracy Relative accuracy.
     * @param absoluteAccuracy Absolute accuracy.
     */
    public BrentSolver(double relativeAccuracy, double absoluteAccuracy) {

        super(relativeAccuracy, absoluteAccuracy);
    }

    /**
     * Construct a solver.
     *
     * @param relativeAccuracy Relative accuracy.
     * @param absoluteAccuracy Absolute accuracy.
     * @param functionValueAccuracy Function value accuracy.
     */
    public BrentSolver(double relativeAccuracy, double absoluteAccuracy,
                       double functionValueAccuracy) {

        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
    }

    /**
     * {@inheritDoc}
     */
    @Override
    protected double doSolve() {

        double       min                   = getMin();
        double       max                   = getMax();
        final double initial               = getStartValue();
        final double functionValueAccuracy = getFunctionValueAccuracy();

        verifySequence(min, initial, max);

        // Return the initial guess if it is good enough.
        double yInitial = computeObjectiveValue(initial);

        if (FastMath.abs(yInitial) <= functionValueAccuracy) {
            return initial;
        }

        // Return the first endpoint if it is good enough.
        double yMin = computeObjectiveValue(min);

        if (FastMath.abs(yMin) <= functionValueAccuracy) {
            return min;
        }

        // Reduce interval if min and initial bracket the root.
        if (yInitial * yMin < 0) {
            return brent(min, initial, yMin, yInitial);
        }

        // Return the second endpoint if it is good enough.
        double yMax = computeObjectiveValue(max);

        if (FastMath.abs(yMax) <= functionValueAccuracy) {
            return max;
        }

        // Reduce interval if initial and max bracket the root.
        if (yInitial * yMax < 0) {
            return brent(initial, max, yInitial, yMax);
        }

        throw new NoBracketingException(min, max, yMin, yMax);
    }


    /**
     * Search for a zero inside the provided interval.
     * This implementation is based on the algorithm described at page 58 of
     * the book
     * <quote>
     * <b>Algorithms for Minimization Without Derivatives</b>
     * <it>Richard P. Brent</it>
     * Dover 0-486-41998-3
     * </quote>
     *
     * @param lo Lower bound of the search interval.
     * @param hi Higher bound of the search interval.
     * @param fLo Function value at the lower bound of the search interval.
     * @param fHi Function value at the higher bound of the search interval.
     * @return the value where the function is zero.
     */
    private double brent(double lo, double hi, double fLo, double fHi) {

        double       a   = lo;
        double       fa  = fLo;
        double       b   = hi;
        double       fb  = fHi;
        double       c   = a;
        double       fc  = fa;
        double       d   = b - a;
        double       e   = d;

        final double t   = getAbsoluteAccuracy();
        final double eps = getRelativeAccuracy();

        while (true) {

            if (FastMath.abs(fc) < FastMath.abs(fb)) {

                a  = b;
                b  = c;
                c  = a;
                fa = fb;
                fb = fc;
                fc = fa;
            }

            final double tol = 2 * eps * FastMath.abs(b) + t;
            final double m   = 0.5 * (c - b);

            if ((FastMath.abs(m) <= tol) || Precision.equals(fb, 0)) {
                return b;
            }

            if ((FastMath.abs(e) < tol) || (FastMath.abs(fa) <= FastMath.abs(fb))) {

                // Force bisection.
                d = m;
                e = d;
            } else {

                double s = fb / fa;
                double p;
                double q;

                // The equality test (a == c) is intentional,
                // it is part of the original Brent's method and
                // it should NOT be replaced by proximity test.
                if (a == c) {

                    // Linear interpolation.
                    p = 2 * m * s;
                    q = 1 - s;
                } else {

                    // Inverse quadratic interpolation.
                    q = fa / fc;

                    final double r = fb / fc;

                    p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
                    q = (q - 1) * (r - 1) * (s - 1);
                }

                if (p > 0) {
                    q = -q;
                } else {
                    p = -p;
                }

                s = e;
                e = d;

                if ((p >= 1.5 * m * q - FastMath.abs(tol * q))
                        || (p >= FastMath.abs(0.5 * s * q))) {

                    // Inverse quadratic interpolation gives a value
                    // in the wrong direction, or progress is slow.
                    // Fall back to bisection.
                    d = m;
                    e = d;
                } else {
                    d = p / q;
                }
            }

            a  = b;
            fa = fb;

            if (FastMath.abs(d) > tol) {
                b += d;
            } else if (m > 0) {
                b += tol;
            } else {
                b -= tol;
            }

            fb = computeObjectiveValue(b);

            if (((fb > 0) && (fc > 0)) || ((fb <= 0) && (fc <= 0))) {

                c  = a;
                fc = fa;
                d  = b - a;
                e  = d;
            }
        }
    }
}
